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Abstract. We present a novel interpretation of IceCube high energy neutrino events (with 
energy larger than 60 TeV) in terms of an extraterrestrial flux due to two different contribu¬ 
tions: a flux originated by known astrophysical sources and dominating IceCube observations 
up to few hundreds TeV, and a new flux component where the most energetic neutrinos come 
from the leptophilic three-body decays of dark matter particles with a mass of few PeV. Dif¬ 
ferently from other approaches, we provide two examples of elementary particle models that 
do not require extremely tiny coupling constants. We hnd the compatibility of the theoretical 
predictions with the IceCube results when the astrophysical flux has a cutoff of the order 
of 100 TeV (broken power law). In this case the most energetic part of the spectrum (PeV 
neutrinos) is due to an extra component such as the decay of a very massive dark matter 
component. Due to the low statistics at our disposal we have considered for simplicity the 
equivalence between deposited and neutrino energy, however such approximation does not 
affect dramatically the qualitative results. Of course, a purely astrophysical origin of the 
neutrino flux (no cutoff in energy below the PeV scale - unbroken power law) is still allowed. 
If future data will conhrm the presence of a sharp cutoff above few PeV this would be in favor 
of a dark matter interpretation. 
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1 Introduction 


After more than 80 years from its first evidence in the Coma galaxy cluster by Fritz Zwicky, the 
nature of Dark Matter (DM) still remains an open question. Elementary particle physics can 
provide interesting schemes where to allocate viable DM candidates, one of the most attractive 
and simple scenarios being the Weakly Interacting Massive Particle (WIMP) scenario with 
a DM mass in the range 0 ( 1 ) GeV - 0(100) TeV [1]^ and interaction rates of the order of 
weak interactions. Even though such schemes naturally emerge in the SUSY extension of 
the electroweak Standard Model (SM), up to now almost all indirect or direct searches have 
not provided any clear evidence [3], and DM observations remain linked to their indirect 
gravitational footprint only. 

Indeed, a large amount of different theoretical frameworks has been proposed in litera¬ 
ture, where the mass of DM candidates is spread over many order of magnitude, from about 
10“^^ GeV up to 10^® GeV, e.g. axions, KeV sterile neutrinos, majorons, or the heaviest wim- 
pzilla (~ 10^^ GeV). Lacking any direct DM detection at LHG experiments, it is likely that 
the only viable way to look for very massive DM candidates would exploit indirect searches 
in astrophysical observations. From this point of view, neutrino telescopes like IceGube (IG) 
provide a chance to observe high energy cosmic ray phenomena induced by such massive DM 
particles, where energetic neutrinos are produced. 

The IceGube Neutrino Observatory experiment [4, 5] is an excellent example of high 
energy neutrino astronomy. This rapidly developing branch of physics is important for differ¬ 
ent reasons [6-8]. In fact, by looking at neutrinos, which are neutral and weakly interacting 
particles, one can trace back the sources even in presence of an intergalactic background 
and magnetic fields. Moreover, the observation of astrophysical neutrinos is of paramount 
importance, both because their presence is a proof that acceleration of hadronic matter is 
involved, and because the features in their energy spectrum may indicate their production in 
a top-down framework from interactions involving heavy non standard particles. 

During three years (2010-2013) IG Gollaboration has observed several neutrino events 
in the TeV-PeV range. While the lower energy ones can be well explained in terms of atmo¬ 
spheric neutrinos, events in the range between 100 TeV and 2 PeV seem to be due to some 

^This upper limit arises from a model independent unitary constraint. However, a typical mass for WIMP 
is ~ 0 { 1 ) TeV [2]. 
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extraterrestrial process. There are many astrophysical possible sources for such neutrinos [9], 
both from pp [10-12] and p'j [13] interactions; the bottom-up scenarios range from extragalac- 
tic Supernova Remnants (SNR) [14] to Active Galactic Nuclei (AGN) [15, 16] and Gamma 
Ray Bursts (GRB) [17], all of these with specihc emission spectra due to the different pro¬ 
duction environment. For example, p'j spectrum is peaked [13] while pp spectrum is flatter 
[12]. According to the IG analysis [5], galactic sources alone can’t explain the excess and 
extragalactic sources are necessary. However, even assuming an extragalactic origin, it is not 
straightforward to ht all the data; for example, pj AGN spectrum gives a good description 
of high energy events but does not satisfactory ht lower energy data [16]. On the other side, 
GRB predicted spectra agree with data modulo the normalization of hux, which, unfortu¬ 
nately, has an upper limit (given by searches for correlation with observed GRB) more than 
one order of magnitude below the observed value [18]. 

More recently it has been proposed that high energy IG events could be related to the 
decay of DM particles [19-30]. Note that the presence of PeV decaying DM component can also 
alleviate some tension among cosmological parameters estimates [31]. Such an interpretation 
would be supported by a lack of events above 2 PeV compatible with the decay of a massive 
particle with a mass of this order of magnitude, and by a distribution of arrival directions 
apparently uncorrelated with the galactic disk. If this interpretation is correct, IceGube will 
provide information about the DM mass range and cross-section with ordinary particles and 
thus, important hints for future DM experiments. 

For an elementary particle with a mass of 0{1) PeV, the maximal annihilation cross- 
section obtained by saturating the unitarity bound yields to a completely negligible signal at 
IceGube. One is then compelled to consider decaying DM instead [19]. 

Alternative connections between DM and IG events have been proposed in literature. 
In particular the boosted DM mechanism [32-35] is based on the idea that a highly energetic 
(boosted) population of DM particle is originated by the decay of more massive and long- 
lived non-thermal relic, dominating the DM distribution. Such boosted particles then interact 
with nucleons of detector via neutral current interactions. Finally, different schemes have been 
proposed, where the bump in the neutrino flux at PeV observed by IceGube is explained as 
the s-channel enhancement of neutrino-quark scattering by a leptoquark with a mass of 0.6 
TeV that couples to the r-flavour and light quarks [36]. 

In the present paper, we propose an interpretation of IceGube PeV events in terms 
of leptophilic three-body decays of a DM particle. The operator responsible for the decay 
results to have quite a natural coupling and is selected out by global flavour symmetries. 
Such symmetries forbid at the same time the two body decays and other dangerous three 
body-decay operators. 

The paper is organized as follow. In Section 2 we give a brief review of IG data, and 
in Section 3 we outline the theoretical framework for DM particle. In Section 4 we discuss 
the neutrino flux produced, and in Section 5 we describe our results. Section 6 contains our 
conclusions. 

2 IC data 

The IceGube Neutrino Observatory has been searching during three years (2010-2013) [5] for 
astrophysical neutrinos in the energy range from 30 TeV to 100 PeV. The observed data have 
been analyzed in terms of their energy spectrum, arrival direction, and flavour. The expected 
background arises from muons and neutrinos coming from the decays of vr and K produced by 
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cosmic ray interactions in the atmosphere. Actnally, when the energy of the parent meson in¬ 
creases, its lifetime becomes longer and correspondingly, the interaction probability dominates 
on the decay, giving a snppression in the atmospheric mnon and nentrino flnx at high energy. 
No snppression is instead foreseen for the atmospheric nentrino backgronnd coming from the 
decay of charmed mesons, the so-called prompt component, becanse of the short lifetime of 
these mesons. Measnring the mnon detection rate in a separate region of the telescope, the 
IC collaboration gives the following estimation of the mnon and nentrino backgronnd 

N^± = 8.4 ± 4.2 , (2.1) 

where stands for the nnmber of all flavonr nentrinos and antinentrinos and its asym¬ 

metric error is dne to the prompt component. 

IceCnbe collected 37 events in 988 days, with deposited energies ranging from 30 TeV 
to 2 PeV. In particnlar, three events were detected, with deposited energy of the order of 
PeV, which are the most energetic nentrino events ever detected. Among all the events, two 
of them, events 28 and 32, having snb-threshold signals in IceTop, seem to be part of the 
expected mnon backgronnd (in particnlar, event 32 cannot be reconstrncted with a single 
direction and energy), while three ambignons downgoing tracks seem to be of an atmospheric 
origin. 

It has been argned in literatnre [37] that, dne to the present nncertainty on the overall 
rate and starting energy of the prompt nentrino component, by assnming a different cross 
section for charmed meson prodnction and/or a slightly different cosmic ray primary flnx one 
conld modify the expected atmospheric nentrino backgronnd in snch a way to rednce the 
necessity of an extra nentrino flnx in the high energy range. In this case since abont 50% of 
downgoing prompt nentrinos arrive together with mnons, which shonld trigger the mnon veto, 
npgoing events shonld be more than downgoing ones at these energies (“sonthern hemisphere” 
snppression). IceCnbe observes exactly the contrary, so it presently seems nnlikely that the 
atmospheric backgronnd and the observed nentrino flnx conld be reconciled thanks to prompt 
nentrinos [38, 39]. 

Another indication comes from the topology of detected events, which belong to two 
classes: track events, associated with the propagation of a high energy mnon, and shower 
events, which correspond to the prodnction of a large aggregate of secondary particles (with 
similar energy resolntion in the two cases bnt, of conrse, better angnlar resolntion in the first 
one). By assnming that data are dne to a pnrely conventional atmospheric flnx, one shonld 
connt more tracks than showers (since atmospheric nentrinos are mainly mnon nentrinos). 
On the other side, a (1/3 : 1/3 : 1/3) flavonr proportion in the nentrino flnx, as expected for 
astrophysical nentrinos, wonld resnlt in only 20% CC interactions, a closer resnlt to the 
IceCnbe finding of less tracks (24%) than showers (76%). 

On the basis of this analysis the IC collaboration conclndes [5] that a pnrely atmospheric 
origin of the high energy events, reqniring a 3.6 times higher charm normalization, is rejected 
at 5.7cJ. Moreover, the data are well described in terms of a global fit inclnding backgronnd 
atmospheric mnons and nentrinos, prompt nentrinos and an isotropic astrophysical flnx, which 
for each flavonr takes the form (qnoted errors are l-cr nncertainties) 

^ 2 ^;^ = (0.95 ± 0.3) X 10"®GeV cm’^ sr"^ . (2.2) 

dE 

This resnlt satisfies the Waxman-Bahcall bonnd for optically thin sonrces [40], obtained snp- 
posing that all the charged particles created by cosmic accelerators give their energy to kaons 
and pions. 
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While the unbroken E~'^ hypothesis gives a fairly good description of the data energy 
spectrum, it predicts 3.1 more events at 2 PeV, seeming to require a softer spectrum or a high 
energy cutoff. A more general model of the astrophysical component by a piecewise function 
of the energy gives the best-fit 

^2 = 1.5 X 10-« ( GeV cm-2 s"! sr'^ . (2.3) 

dE V100 TeVy ^ ^ 

This corresponds to the lower boundary of the total statistical and systematic uncertainty on 
the energy power law (with a zero charm contribution), which on the other side reaches, at 
90% C.L., the previous E~‘^ behaviour. 

In the present work, by following the background analysis of ref. [5], we assume that 
neutrino flux is mainly dominated by the atmospheric component up to 60 TeV, whereas 
we consider a bottom-up neutrino contribution from known astrophysical sources (like for 
example extragalactic SNR) in the [60 TeV, 300 TeV] range (hereafter denoted as astrophysical 
neutrino flux), and, at the same time, a top-down additional component at higher energy 
(hereafter denoted DM neutrino flux). This would naturally produce a sharp cutoff, observed 
in the data around few PeV, and if confirmed by new data, represents an intriguing feature 
of the IC observations. 

3 The Theoretical Framework 

For a heavy fermion singlet DM candidate, Xi which is directly coupled to neutrinos, the 
lowest dimension coupling is a Yukawa interaction 

yLa^X , (3.1) 

where cj) = ia24>* {4’ denoting the SM Higgs doublet), La=e,^i,T stand for the lepton doublets, 
and (T 2 is a Pauli matrix. As shown in [23], this new interaction term is cosmologically 
safe (sufficiently long lived DM), and relevant for IceCube (a DM particle specie abundant 
enough to fit the observed ffux) for a fine tuned tiny coupling, y ~ 0(10“^^). Note that the 
interaction term in Eq. (3.1) is reminiscent of the right-handed neutrino in see-saw models, 
but its contribution to the lightest neutrino mass would be negligible due to the smallness of 
the y coupling. The operator in Eq. (3.1) would yield a sharp peak in energy. The presence of 
the Higgs held in Eq. (3.1) allows for an abundant production of secondary neutrinos via the 
decay of the Higgs particles to heavy quarks, giving an almost hat neutrino hux at energies 
lower than PeV [23]. Such a result could turn out to be problematic if known astrophysical 
contributions were included in the analysis. For instance in Ref. [14], low energy (up to 100 
TeV) IC events are shown to be well htted by means of standard extragalactic SNR. 

In the present paper we try both to improve the need of an unnatural coupling and at 
the same time, to reproduce the IC data including an astrophysical neutrino component. This 
suggests to consider an interaction term that does not directly involve quarks, the Higgs held 
or gauge bosons (leptophilic DM), and which is a higher dimension operator. In this case, 
the feeble coupling can be understood in terms of a large mass scale. In particular we assume 
that: 

1) the DM held x is coupled to the SM particles via a leptophilic coupling; 

2) the lifetime of x is suppressed by powers of the scale of new physics entering in the 
non-renormalizable coupling; 
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Dimensions 

DM decay operators 

4 


5 

— 

6 

LI Lx, ^^Ux, 


QdLx, uQLx, LdQx, u^f^dh^x, 


D>^^D^Lx, D^^D^^Lx, 




Table 1. Gauge-invariant operators up to dimension-6 inducing fermion singlet DM decay as taken 
from Table 1 of Ref. [41]. The adopted notation is explained in the text. 


3) there is a direct coupling to neutrinos allowing for a primary neutrino flux with energy 
of the order of y mass. The multi-body final state may contain lower energy neutrinos 
as well, so that neutrino flux also spreads to lower energies; 

4) mass and couplings of x are determined in order to produce a neutrino flux just domi¬ 
nating the energy region around PeV. This implies, on the contrary, that the flux in the 
energy region [60 TeV, 300 TeV] is due to astrophysical sources as will be discussed in the 
following. We also assume for simplicity that x represents the dominating contribution 
of Cold DM. 

As in [41], one can list the gauge-invariant operators up to dimension- 6 , shown in Table 1. L 
and Q are the lepton and quark weak doublets and u and d the corresponding right-handed 
quarks. The field i stands for the right-handed lepton singlet, and finally, Wfj,^ and are 
the field strength tensors of SU{2)l and U{1)y gauge bosons. To simplify notation, we have 
omitted the family index for matter fields. 

Remarkably, the only operator in this list which satisfies requirements l)-3) is the non- 
renormalizable lepton portal: 


|aa {L,x) + h.c.. (3.2) 

IVlpi 

where {a, /3, 7 } are flavour indices thus labeling 27 operators, and the round brackets indicate 
the Lorentz contractions. The mass scale here is chosen as the Planck mass Mpp Indeed, for 
this choice the order of magnitude of couplings y will not be unnaturally small, see Section 5. 

If this operator is the only source of DM decay, one has to invoke some selection rule 
which forbids the dimension-4 operator, which, as we mentioned, is compatible with IC results 
for an extremely small coupling only. This can be done by using global flavour symmetries, 
both Abelian like Uf{l) and non-Abelian groups like A 4 , A(27), etc. We will focus in the 
following on two benchmark schemes (see Appendix for more details about these models)^: 

• model 1 ) - symmetry {a, /3,j} = {fi,e,T} {T,e, y} + {e, 

^In our analysis we assume for simplicity that y is a Dirac fermion. If x was a Majorana fermion, one 
could not reproduce, for example, the Abelian model 1) that requires a flavour charge for x different from 0, 
whereas a scheme like A 4 would still be possible. 
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• model 2) - A 4 symmetry {a,/3,7} = {e,/i, r} + cyclic permutations; 

where the brackets show the flavour assignments corresponding to non vanishing couplings 
2 / 0 / 37 - Note that expanding the SU{2) contractions, the operator of Eq. (3.2) always yields a 
coupling of DM with two charged leptons and one neutrino. Depending on the flavour index 
the charged leptons can then possibly decay producing secondary neutrinos. 


4 DM neutrino flux 


Following the method outlined in Ref. [23] we consider both the contributions to top-down 
neutrino flux coming from the galactic and extragalactic distributions of DM particles x (the 
sum of X and x contributions is implicitly assumed) 






djG dJ/^G 


dE„ 


dE„ 


(4.1) 


where the solid angle integration is on the longitude and latitude in the galactic coordinate 
system, I and b. In the following we will always sum the flux of neutrinos and antineutrinos 
of different flavours, unless explicitly stated. 

The galactic component due to the Milky Way halo can be written as 


d 


1 


47r My. r- 


X 'X 


cx=e^fi,r 


dN^+, 

dE^ 


rc 

(E.) / 

Jo 


ds p^{r{s,l,b)) 


(4.2) 


where and denote the mass and lifetime of DM particle. The quantity Py.{r) denotes 
the density prohle of DM particles in our Galaxy as a function of distance from the Galactic 
center, r, and dN^_^_p/dEi, stands for the energy spectrum of neutrinos and antineutrinos of 
flavour a produced in the decay. Note that the parameter s is related to r via the expression 


r{s,l,b) = \J+ Rq — 2sRq cos b cos I 


(4.3) 


where Rq ~ 8.5 kpc is the distance of the Sun from the Galactic center. In the following we 
assume a Navarro-Frenk-White density prohle for the halo 


Pxir) = 


Px 


r/rc(l -h r/vcY 


(4.4) 


with Tc = 20 kpc the critical radius and py. = 0.33 GeV cm“^. We have checked the dependence 
of our predictions on the choice for the density prohle (Einasto, Isothermal, etc.). The results 
are very similar, with a variation of the best ht values of the order of few percents. 

The quantity dN^^p/dEy has been evaluated by means of a MonteGarlo procedure. 
In our code we take into account primary neutrinos (antineutrinos) and also the neutrinos 
(antineutrinos) produced by the decay of p and r leptons. The r leptons have different decay 
channels that involve pions as well [3]. The MonteGarlo takes into account ~ 90% of the r 
decay width, including the two leptonic decays into muons and electrons, and semileptonic 
decay channels up to three pions in the hnal state, whose charged states eventually decay 
mainly producing p and corresponding neutrinos (antineutrinos). It is worth observing that 
the electroweak radiative corrections produce a large number of secondary neutrinos that are 
typically placed at energies almost two orders of magnitude smaller than the mass of DM 
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particle [42], Hence, for a mass of few PeV, such low energy neutrinos do not significantly 
affect our analysis that is restricted to neutrino energies larger than 60 TeV. In our study, 
we also consider the total energy injected by DM decay in the electromagnetic (e.m.) sector, 
and compare such a value with the bound provided by Fermi-LAT [43]. 

The isotropic extragalactic component of the differential flux is given by 


dE^ 






dz 


E 


div; 




H{z) ^ 

'■ ' a=e,fi,T 


dE, 


((1 + ^w 


(4.5) 


where H{z) = + Hm(l + z)^ is the Hubble expansion rate as a function of redshift 

z and Per = 5.5 x 10“® GeV cm“^ is the critical density of the Universe. We assume a 
ACDM cosmology with parameters Ha = 0.6825, Hm = 0.3175, = 0.2685 and h = 

i7o/100kms“^ Mpc“^ = 0.6711, according to Planck experiment results [44]. We found that 
the galactic and extragalactic components of the DM differential flux are of the same order 
of magnitude. 

The x-lifetime, r^^., depends on the values assumed for yap'y- For models 1) and 2) we 
get 

model 1) ^3 \^\y^leT - 2/re/i| + \ye^ie\ j , (4.6) 

model 2) (-1-) 

where and y_ denote the two independent symmetric and antisymmetric couplings in A 4 
[41]. In the previous expressions we neglected all final state masses, due to the large value of 
M^. 

In the following, for simplicity we assume real couplings and the relations 


model 1) - We/dl = \yefie\ = 2 / , (4.8) 

model 2) |y+| = |y_| = y . (4.9) 


4.1 DM relic abundance 

As well known, the standard thermal freeze-out mechanism is not a viable option for particles 
with masses larger than 0(100) TeV because of the unitarity bound on the cross-section [1]. 
Hence, considering a DM candidate of PeV mass raises the question of the production mecha¬ 
nism, which is crucial to determine its relic abundance. There exist several scenarios, beyond 
the WIMP paradigm, which make such a heavy particle a perfectly viable DM candidate (see, 
e.g.. Ref. [45] for a review). 

One possibility is to consider a non-thermal production during reheating of the Universe, 
after the inflationary stage, in a cosmological scenario with a low reheating temperature [46, 
47]. After inflation, the Universe reaches its maximal temperature, Tmax- This temperature 
can be much higher than the reheating temperature, Tr//, defined as 

X 1/4 

] a/TVWh , (4.10) 

where is the effective number of relativistic degrees of freedom at Trr, and 

Pp the inflaton held decay rate. If the y particle production takes place between Tmax and 


Trh = 0.2 


200 
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Trh, the DM abundance scales as [46]: 



■3/2 



~ M^{av) 


(4.11) 


where (av) is the thermal average of the DM annihilation cross-section times the M0ller flux 
factor. Notice that the x particles are never in chemical equilibrium (even if they are in kinetic 


equilibrium), so their abundance is described by a power law instead of being exponentially 
suppressed as exp(—Assuming (av) ~ obtain for = 200 and = 


IPeV 


~ 370GeV , 


(4.12) 


if X provides the whole contribution to DM today. The previous equations are valid under 
the condition Tmax > > 20 Tjih [48], corresponding to the fact that DM never reached 

local thermal equilibrium in the early Universe (which would otherwise spoil the result of 
eq. (4.11)). From eq. (4.12) we see that the lower bound is readily satisfied. 

Another viable possibility is the production through inelastic scattering interactions 
among the thermal plasma and energetic particles originating from inflaton decay and taking 
place at temperatures between Tmax and Tjih [49]. In this case, assuming that DM couplings 
are of same order as gauge couplings and that Mp > M^/2 Tjih, with Mp the mass of the 
inflaton, in order to account for the observed DM abundance the following condition should 
be satisfied [49]: 



(4.13) 


This mechanism would require quite low, thought viable, values for the reheating temperature 
to account for ~ 1 PeV. 

Finally, PeV dark matter could also be produced from inflaton decays (directly or through 
cascades). However, this mechanism is quite model-dependent because it strongly depends on 


the couplings of the inflaton to the particles of the model [50, 51]. For the sake of illustration, 
the direct production of DM from inflaton decay gives the abundance (assuming all decays 
happen at Trh): 


^ 2.1 X 10* DM) , (4.14) 


(4.14) 


where B{p ^ DM) is the average number of DM particles produced per inflaton decay. 

5 Analysis and results 


In order to reproduce the sharp cutoff observed in the IC neutrino data, in the present 


analysis we assume that the neutrino flux is given by the combination of an atmospheric 
component up to 60 TeV [5] and two different neutrino contributions, a bottom-up one from 
known astrophysical sources, like extragalactic SNR, in the [60 TeV, 300 TeV] range (hereafter 
denoted as JAst); and a top-down component at higher energies (that is the flux J^, defined 
in Eq. (4.1)). Hence, for > 60 TeV the total neutrino flux is given by 
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Then, the number of neutrinos in a given energy bin [Ei, -Ej+i] is equal to 


Ni = 47rAt 



Aa (E) , 


(5.2) 


where At = 988 days is the exposure time after 3-years of IC experiment, and Aa (E) is the 
neutrino effective area [4] for different neutrino flavour a. In order to compare the theoretical 
predictions of Eq. (5.2) with the observations we need to recast such expression in terms of 
deposited energy. In general, to statistically estimate the ratio between the deposited and 
neutrino energies a MonteCarlo simulation of the apparatus is required [52]. When for a bin of 
the deposited energy a significant statistics is collected one could apply an average ratio that 
results to be of the order of 97%-|-23 %)/(cj‘"‘"-|- o"^^) ~ 75% (see Table 1 of [52]). 
Remarkably, such a number appears to be quite stable as a function of the neutrino energy. 
Unfortunately, due to low statistics collected till now, this procedure would be characterized 
by a large uncertainty in the energy determination. For this reason we prefer in this paper 
to assume the simplicity ansatz that the two energies coincide. Notice that in any case an 
expected shift in the energy of the order of 25% is not going to change dramatically our 
conclusions. 

Known astrophysical sources able to produce a neutrino flux similar to the one observed 
by IceCube are high energy accelerators [9] such as extragalactic SNR [14], Hypernova Rem¬ 
nants (HNR) [14], AGN [15, 16], and GRB [17]. It is worth observing that, due to the large 
uncertainties in the parameters related to the physics of the accelerator mechanisms, there is 
enough room to reproduce IG neutrino flux. In particular, an extragalactic SNR can produce 
a neutrino flux with a cutoff 0(100) TeV coming from the proton-proton hadronic interactions 
up to few PeV. On the other hand the HNR are sources of 100 PeV protons that can explain 
the PeV neutrino events [14]. While classical GRB are ruled out, low-luminosity GRB and 
choked jets, which cannot be triggered by GRB satellites, are totally viable as the origin of 
PeV neutrinos [53-55]. AGN can be the sources of PeV neutrinos [16, 56], but if the flux is 
normalized to match the observed IG PeV events, the lower energy part is inconsistent with 
data. All these astrophysical sources have also an associated 7 -ray flux that would give a 
significant contribution to the diffuse 7 -ray background, strongly constrained by Fermi-LAT 
[43]. 

In order to parametrize the astrophysical flux one can use either an Unbroken Power 
Law (UPL), with a power law behavior in the whole IG region, or a Broken Power Law 
(BPL) where an exponential cutoff is assumed at some energy scale Eq. Gonsidering both 
options essentially covers the wide range of accelerator mechanisms related to the different 
astrophysical sources. Hence, using the notation adopted by the IG Gollaboration, we will 
consider: 


i) Unbroken Power Law (UPL): 

2 d JAst 

ii) Broken Power Law (BPL): 

2 dJAst 


{Eu) = Jo 


E„ 


100 TeV 
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(Eu) = Jo 


E, 


100 TeV 


exp 


E, 

Eo 


(5.3) 


(5.4) 
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Model 

Case 

y [10-"] 

Jo [10-«] 

xVdof 

1) Uf{l) 

UPL 

BPL 

1 0+°-^ 
-‘-■^-0.7 

1 1+0.6 
-‘-•-‘^-0.5 

0 t::'+2.8 

10.3/12 

9.2/12 

2) A4 

UPL 

BPL 

f) Qr;+0.21 

U.OO_o 21 

Q oy+0.17 
-0.16 

2.4t|8 

10.7/12 

9.6/12 


Table 2. The marginalized 95% C.L. for parameters y and Jq (expressed in unit of GeV cm“^ s“^ 
sr“^) corresponding to models 1) and 2), and UPL and BPL parameterizations, respectively. The last 
column reports the reduced 


where 7 + 2 is the spectral index and Jq is the flux normalization. In the present analysis we fix 
the value of Eq to be equal to 125 TeV in agreement with the extragalactic SNR results [14]. 
Furthermore, we restricted the spectral index to the physical range 7 G [0,1], as suggested by 
cosmic accelerator mechanisms. 

The fit has been done by means of a multi-Poisson likelihood analysis [58] in which the 
takes the expression 
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where Ni is the expected number of neutrinos for energy bin provided by Eq. (5.2), while n* is 
the observed one. Once the atmospheric background has been subtracted, we fit the IC data 
by using the parametrization of the astrophysical flux in both cases of UPL and BPL, and 
considering model 1) and model 2) of Section 3 for the DM neutrino flux. The mass of DM 
particle has been varied in a range [1 PeV,10 PeV] able to produce a drop in the flux. The 
best fit value is found for = 5.0 PeV independently of the model adopted. On the other 
side, after scanning the 7 range, we find a best fit value 7 = 1.0 for UPL (spectral index = 
3.0), and 7 = 0.0 for BPL (spectral index = 2.0). 

In Table 2, we give the marginalized 95% C.L. for parameters y and Jq, for models 1) 
and 2), and UPL and BPL parameterizations, respectively. From the values of reduced 
shown in the last column, we see that the experimental data slightly prefer the BPL scheme 
with respect to the UPL one. In each case, the analysis shows that a non-vanishing DM 
contribution at 2-a level (95% C.L. for y not compatible with zero) is required. This is 
mainly due to the presence of the sharp cutoff in the data at high energy. By comparing the 
results obtained for Uf{l) and +.4, one cannot appreciate a significative difference. Indeed, 
the two models essentially provide similar features at the level of produced neutrino flux. For 
the cases reported in Table 2 we have checked that the total energy injected by DM decay in 
the e.m. sector is smaller than the bound provided by Fermi-LAT [43]. The neutrino events 
as function of E^, for model 1) are shown in Fig. 1 (first row) and compared with IC data 
for the two models, DM + UPL (column A) and DM+BPL (column B). In the second row 
we report the 68 % C.L. (dashed) and 95% C.L. (solid) contours for the two parameters, y 
(defined in Eq. (4.8)) and Jq (defined in Eq.s (5.3) (5.4)), in the cases DM + UPL (column 
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Figure 1. Results of the analysis for model 1). First row shows the neutrino events as a function 
of the neutrino energy Ei, for the DM + UPL (column A) and DM+BPL (column B) models. The 
red (long-dashed) line is the best fit (background + astrophysical + DM components), and its band 
represents the 68 % C.L. resulting from the fit. The purple (dashed) and green (solid) lines are the 
astrophysical and DM contributions, respectively. The black points are the IC data, and the blue 
region shows the upper limit for the sum of all backgrounds (see Ref. [5]). In the second row we report 
the 68 % C.L. (dashed) and 95% C.L. (solid) contours for the two parameters y and Jq corresponding 
to DM + UPL (column A) and DM+BPL (column B). The crosses are the best-fit points. 


A) and DM+BPL (column B). The crosses stand for the best-fit points. The same quantities, 
but corresponding to the model, are reported in Fig. 2. The plots in Figs. 1 and 2 show 
that there is no significative difference between the two models: both two models predict 
the observation of neutrinos in the energy range [0.3 PeV, 1.0 PeV] and a sharp cutoff at 
the energy of few PeV. Notice that since in both (7/(1) and symmetry cases the galactic 
and extragalactic components of the DM neutrino flux are of the same order of magnitude, 
we expect at the PeV energies an almost isotropic neutrino flux with a significant level of 
anisotropy near the galactic center. This is in a good agreement with the IC observations. 

In order to understand the reason for the strong similarity between the predictions of 
model 1) and 2), we have also considered two situations where the operator of Eq. (3.2) is 
characterized by a single term only, {a, /3, 7 } = {e, e, e} and {a, /3, 7 } = {r, r, r}, though 
they cannot be realized by using the flavour symmetries considered here, as shown in the 
Appendix. From a phenomenological point of view it is interesting to consider such cases 
because they correspond to the extreme situations with respect to the number of r leptons 
produced and thus, of secondary neutrinos from their decays. In the left panel of Fig. 3 we 
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Figure 2. Same as Fig. 1 for model 2). 




Figure 3. Left panel shows the flavour compositions at Earth of the DM neutrino flux for model 1), 
2) and the two fully diagonal cases (e, e, e) and (r, r, r), as well (see text). The green square represents 
the IC flavour analysis of Ref. [57] , but referred to the integrated neutrino flux above 35 TeV. The 
prediction for models 1), and 2) and the (r, r, t) case are represented by the red disk, whereas the 
blue star stands for the (e, e, e) case. Right panel presents the DM neutrino flux for model 1) (green, 
solid), model 2) (purple, solid), (e,e, e) (blue, dashed), and (r, r, t) (red, long-dashed). The two 68 % 
C.L. bands refer to the two cases (e,e, e) and (r, r, r). 
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report the flavour compositions at Earth of the DM neutrino flux for model 1) and 2) and the 
two fully diagonal cases (e, e, e) and (r, r, r), as well. The black dot represents the IC flavour 
analysis of Ref. [57]. Models 1) and 2) and the (r, r, r) diagonal case provide the same flavour 
composition at Earth, namely (/e : : fr)^ ~ (1/3 : 1/3 : 1/3), which are the flavour ratios 

provided by the standard astrophysical sources. A different situation is given by the diagonal 
case (e, e, e) where we have only electron (anti)neutrinos as DM decay products, that leads at 
Earth (/e : : fr)^ ~ (0.55 : 0.19 : 0.26). It is worth reminding that a study of shower/track 

composition of the PeV events can shed light on the flavour ratios and hence on the flavour 
structure of the coupling. Such analysis cannot be performed at the moment due to the low 
statistics available. In the right panel of Fig. 3 one can see that changing the flavour structure 
of the DM-SM coupling does not appreciably affect the DM neutrino flux. 


6 Conclusions 

During the last three years of observation the IceCube Collaboration has observed neutrino 
events in the TeV-PeV range. Even though lower energy events can be well explained in terms 
of atmospheric neutrinos, the most energetic processes detected, in the range between 100 TeV 
and 2 PeV, seem to originate from some extraterrestrial source. In the present analysis, by 
following the results on the background presented in Ref. [5], we assumed an extraterrestrial 
neutrino flux dominating for energies larger than 60 TeV. In particular, in order to recover the 
drop of the flux observed around few PeV, this high energy flux is obtained as the sum of two 
different components: a bottom-up neutrino contribution coming from known astrophysical 
sources (like for example extragalactic SNR) in the energy region [60 TeV, 300 TeV], and a 
top-down additional component, dominating at higher energies, originated from the decay of 
DM particles with mass of few PeV. The astrophysical flux was parametrized using either an 
unbroken power law or a broken power law with an exponential cutoff. We have considered 
both the options UPL and BPL in order to cover the wide range of possible astrophysical 
sources. The top-down term in the neutrino flux is produced by the decay of a PeV DM 
particle, X: which we assume for simplicity to be the dominant cold DM term. The decay of 
y has been calculated by means of a MonteCarlo procedure. 

A similar approach was already investigated in literature [23], but assuming the presence 
of a trilinear coupling for x like the one reported in Eq. (3.1). In this case, the y-decay 
induces a shower whose hadronic content yields secondary neutrinos that provide an almost 
flat behavior of neutrino spectrum at low energy. The main problem of such an approach 
is that it could be in contrast with the standard astrophysical sources since the trilinear 
interaction term alone totally explains the IC extraterrestrial neutrino flux. Moreover, in 
this scenario an unnatural tiny coupling for the interaction term of Eq. (3.1), 0(10“^°), is 
required. In the present analysis we use the same approach of [23], but assume the presence of 
special flavour symmetries (Abelian and non-Abelian), like Uf{l) and A 4 for example. Such 
symmetries forbid a trilinear interaction term of the type in Eq. (3.1), and allow quadrilinear 
leptophilic interaction terms only, like the one in Eq. (3.2). At the same time they reproduce 
the values of neutrino masses and mixing. The leptophilic characteristic is crucial in forbidding 
the production of a huge hadronic component in the decay shower that would flatten again 
the neutrino spectrum. Such a model allows one to improve the problem of naturalness of 
the coupling, since the non-renormalizable terms is weighted by the square of the new physics 
mass scale, that can be as large as the Planck mass. In this scenario the coupling required to 
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recover the IC data results to be 25 orders of magnitude larger than the one for the trilinear 
term. 

The comparison of the prediction for the neutrino events in UPL + DM and BPL + DM 
with the experimental observations allows the determination via a multi-Poisson likelihood 
approach of the free parameters of the model. In Table 2 we show the resulting parameters 
for the Uf{l) and A 4 symmetry models. Both schemes provide a fair description of the 
experimental situation within the statistical uncertainty, with the BPL parameterizations 
providing a slightly smaller value of the reduced The contour plots seem to suggest at 
2 -cj the presence of a hard component (top-down term) in the spectrum (y not compatible 
with zero), whereas the astrophysical contribution could be even vanishing (Jq compatible 
with zero). This is almost expected because the astrophysical contribution partially overlap 
with the atmospheric background (affected by a large uncertainty) up to few hundreds TeV. 
Nevertheless, we cannot conclude at this level of statistics that the presence of DM is really 
demanding. A confirmation by future data of the cutoff above few PeV is therefore particularly 
important. Moreover, since the galactic and extragalactic components of the DM neutrino flux 
are found to be of the same order of magnitude, another evidence in favor of a DM component 
would be the presence of a significant level of anisotropy near the Galactic Center. 

We also analyzed the flavour ratios of neutrinos at Earth predicted by the two models, 
and found that for both models, (/e : : fr)^ ~ (1/3: 1/3: 1/3). Only a coupling of x 

with first lepton generation only would produce clear features in the flavour ratios, which in 
principle could be measurable if a large statistics of events were available. 

Finally, we would like to highlight that our analysis regards only the three years IC data. 
The recent new IC data provide a track event with a deposited energy of about 2.6 PeV, that 
could be explained by increasing the mass of DM particle in our model. However, since such 
an event is not fully contained, the energy of the primary neutrino is unknown and could also 
be very high (0(10) PeV). Note that if the energy of the primary neutrino was of the order 
of few PeV this would not change our conclusions, differently a neutrino energy of the order 
of 10 PeV would be in tension with our analysis. 
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A Flavour symmetry schemes 

As discussed in Section 3, we look for symmetry motivations able to single out the operator 
of Eq. (3.2) among those listed in Table 1. The flavour symmetry schemes that can be used 
may be Abelian or non-Abelian. We will discuss in the following some possible benchmark 
models as relevant examples. 

B Abelian case 

Let us denote with the Uf{l) flavour charge of a generic field i/- R can be shown that 
a flavour Abelian symmetry cannot single out a flavour-diagonal operator in Eq. (3.2). In 
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Table 3. Charges, q^, for the t//(l) symmetry (upper Table). In the lower Table we report the Irre¬ 
ducible Representations allocating SM fields and the dark matter x in the non-Abelian A 4 symmetric 
model (quarks are all singlets under A 4 ). 


fact, in this case both the operators L^iaLaX La<j)(.a should be invariant under such a 
symmetry, but this would imply the invariance of the term La4>X well. Indeed, in terms 
of abelian charges we would have — 2 qL^ + Qx ~ ^ ~QLa +?</> + Qia ~ which 

implies —qLa ~ Qx ~ However, if we mix different lepton flavours then we have more 
freedom to consistently define the charges of a Uf{l) flavour symmetry. This can be done in 
such a way that only a leptophilic dimension 6 operator of Table 1 results to be invariant. A 
possible realization of such a scheme is shown in Table 3. In this case the only invariant DM 
decay operators are 

— i\/f2 ^fJ-^e^rX T Vreiil-'T^eLfiX T Vefie T h.C. . (6.1) 

Mpi 

The charge assignment in Table 3 is by no means unique; different choices would select 
different operators. Moreover with the charge assignment given in Table 3, the dimension five 
Weinberg operator Lo,L^4>(p has a structure very similar to the so-called B 4 two-zeros texture 
given in [59] that can fit the lepton mixing parameters, see for instance [60, 61]. 

In principle the form of the operator is dictated by experiment. A hypothetical study of 
the flavour composition of PeV neutrinos in IC data would give useful hints in the definition 
of the possible flavour charges. 

C Non-Abelian case 

Another possible realization of our scheme invokes non-Abelian discrete symmetries. A simple 
model would employ for instance A(27) by assigning L, (p and £ to 3 Irreducible Representation 
while X remains blind to the symmetry. However, for the sake of definiteness we will consider 
in the rest of this article a framework based on A 4 flavour symmetry with L and i transforming 
as 3 while all the remaining fields are singlets (for more details about such a model we refer to 
[41]). Our conclusions are not affected by this choice. The Lagrangian of the model contains 
the following relevant terms: 

-^(M^FX+h.c.) + 06 , (6.2) 

where Cdim=Q is the lowest order non-renormalizable operator allowed by the symmetries of 
the model. It represents a Fermi-like decay interactions give by 

^+Mp; ^ 

^Note however that this conclusion could be bypassed by Invoking supersymmetry. 
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where the notation (■•) 3 ( 3 ') reflects the fact that in A 4 there are two possible contractions 

of two triplets into one triplet representation. We see that Eq. (6.3) not only reprodnces 

Eq. (3.2) bnt it also signihcantly simplihes its strnctnre. Note that in this case we have two 

independent conplings, namely y and y'. 
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